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ABSTRACT 

We perform a joint fit to differential number counts from Spitzer's MIPS and Her- 
scheFs SPIRE instruments, and angular power spectra of cosmic infrared background (CIB) 
anisotropics from SPIRE, Planck, the Atacama Cosmology Telescope, and the South Pole 
Telescope, which together span 220 < u / GHz < 4300 (70 < A//im < 1400). We simul- 
taneously constrain the dust luminosity function, thermal dust spectral energy distribution 
(SED) and clustering properties of CIB sources, and the evolution of these quantities over 
cosmic time. We find that the data strongly require redshift evolution in the thermal dust SED. 
In our adopted parametrization, this evolution takes the form of an increase in graybody dust 
temperature at high redshift, but it may also be related to a temperature - dust luminosity cor- 
relation or evolution in dust opacity. The counts and spectra together constrain the evolution 
of the thermal dust luminosity function up to z ^ 2.5 — 3, complementing approaches rely- 
ing on rest-frame mid-infrared observations of the rarest bright objects. We are able to fit the 
power spectra without requiring a complex halo model approach, and show that neglecting 
scale-dependent halo bias may be impairing analyses that do use this framework. Our model 
has considerable predictive power and can be used to calculate any one- or two-point statistic 
of the CIB over a wide range of frequency and angular scale. 

Key words: infrared: diffuse background - infrared: galaxies - submillimetre: diffuse back- 
ground - submillimetre: galaxies 



1 INTRODUCTION 

A galaxy's star formation rate (SFR) and f ar-infrared (FIR) emis- 
sion are known to be strongly connected jKennicut^ 1 1998h , due 
to absorption and thermal re-emission of starlight by dusty stellar 
birth clouds. A si mple comparison of the energy in cosmic infrared 
background (CIB: IPuget et al.lll996h photons with infrared emis- 
sion from nearby galaxies indicates that FIR e mission (and thus 
SFR) was considerably higher in the past (e.g. iHauser & Dwekl 
I2OOII and references therein), and extensive measurements of in- 
frared (IR) emission have been made in order to exploit this connec- 
tion and understand how SFR and stellar formation environments 
have evolved over cosmic time. A prevailing picture is that lumi- 
nous and ultra luminous infrared galaxies (LIRGs and ULIRGs; 
10" < Lir/Lq < 10^^ and 10^^ < Lir/Lq < lO", respec- 
tively), which are rare in the local universe, beco me far more im- 



portant to the global SFR at higher redshift (e.g. iLe Floc'h et al 



2005 



2007 



Perez-Gonzalez et al.ll2005l: iDaddi et al.ll2005l : ICaputi et al 



Magnelli et alj|201 ll : iLapi et al.ll201 ih . This trend has pre- 



dominantly been probed using IR luminosity functions (LFs) mea- 
sured using instruments including the Infrared Astronomical Satel- 
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lite (IRAS:lNeugebaueretaLlll984l) . the Spitzer space telescope 



dWerner et al 



2004) an d 



servatorv dPilbratt et allbOlOl) . 



mo re recently, the Herschel Space Ob- 



Current constraints on dust-enshrouded star formation be- 
yond z 2 — 2.5 are limited, for several reasons. Inferring 
the IR luminosity associated with star formation from rest-frame 
mid-infrared measurements (e.g. using Spitzer), is subject to in- 
creasingly large uncertainties at high redshift, relating to, for in- 
stance, lack of information about high-redshift source spectral en- 
ergy distributions (SEDs), and the extent to which high-redshift 
IR emission is associated with obscured active galactic nuclei 
(AGN; e.g.iAlex ander et al. 2005; Lutz et al. 2005; Yan et al. 2003; 
iLe Floc'h etal]l2007l ; ISaiina et al.ll2012l) . Probing the thermal dust 
emission associated with star formation near its rest-frame peak at 
A ~ 100 pm at these high redshifts, using observations at lower 
frequencies, is also difficult. Source confusion, arising from the 
high number density of CIB sources p er instrumental beam area 
(e.g. lBlain eral]|l998l ; lDole et al.ll2003h . greatly impairs efforts to 
resolve the CIB in the submillimeter (submm). Resolved objects 
in maps at 250, 350 and 500 ^^m fro m Herschel's Spe ctral and 
Photometric Imaging Receiver (SPIRE; iGriffin etal.l20 10) account 
for less than a fifth of the total CIB intensity jOliver et al. 201oi). 
Deeper constraints have been obtained with fits to the pixel flux his- 
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togram of confused maps (the 'P{D)' approach: IPatanchon et al 
2009 ; Glenn et al. 2010.), and stacking analyse s (e.g. Dole et al 



IPascale et all l2009l : iBethermin etal 



201 1120 12ch ~ however these studies are typically limited to small 
areas of the sky, and robustly quantifying uncertainties in the 
deep counts is challenging. While submm-selected samples ob- 
served using, for ex ample, the Shared C ommon-User Bolomet- 
ric Array (SCUBA; [Holland et alj 1 19991) have give n us impor- 
tant information about high-redshift dust em i ssion (e.g.lSmail et alj 
19971; iHughes etalj|l998| - lEales et al.|[l999l ; IChapman et alj|2005l ; 



Coppin et alj|2006l) . our understanding is fundamentally limited if 



we have access to only the rare few brightest objects of the popula- 
tion. 

Given the success of recent cosmic microwave back- 
grou nd background (CMB) temperature anisotropy measurements 
(e.g. ISpergel et al.l |2003| ; iLueker et all |2010| ; iFowler et al.l |2010| ; 
iKomatsu et al. 201 ih . it is perhaps not surprising that the an- 
gular power spectrum has emerged as a statistic with which to 
study the propertie s of the numerous unres o lved dusty sources 
(e.g. iLagache et alj I2OO7I; Iviero et al.l |2009|; iHaiian et alj I2OI2I; 



Amblard et alj|201 iLlpTanck Collaboration et alj|201 ll ; IViero et alj 
2012bh . especially given that extragalactic dust emission is a sig- 



nifica nt CMB temper a ture foreground on small angular scales 
(e.g. 'Hall et al j I2OI0I ; iFowler et alj I2OI0I ; iDunklev et all 1201 ll ; 
IShirokoff et alj|201 ID . The fluctuations in the CIB surface bright- 
ness are correlated across angular scales considerably larger than 
a beam, making source confusion largely irrelevant. Dusty galax- 
ies, as with many other luminous populations, are understood to 
exhibit clustering behaviour because galaxies trace the matter den- 
sity field, who se own fluctuations are strongly scale dependent 
jPeebleslflgSol) . Combining an understanding of the dark matter 
clustering with a prescription for how the galaxies and dark mat- 
ter are related, either through a simple biasing prescr iption, or a 



more complex framework, such as the halo mo d el (e.g.lBond et al 



1991b 



20011 ; 



.SeliakI I2OO0I; IPeacock & SmitlJ [2OO0I; IScoccimarro et al, 
Coorav & Shetril2002l) . allows us to constrain the environ- 



mental properties of galaxies (such as characteristic host halo mass 
scales) using measurements of the clustering. The galaxy correla- 
tion function is the statistic that has been widely for these purposes 
in analyses of, for example, local galaxies in the main Sloan Digita l 
Sky Survey (SDSS) sample ( e.g. IZehavi et alj|2002l. |200^. I2OI ll) . 
massive galaxies at z ~ 0.5 l|W hite et al.' 201 1) and luminous red 
galaxies (LRGs; e.g. iBlake et al . 2007; Zheng et al. 20(^. Inter- 
preting angular power spectra of unresolved CIB sources is more 
challenging, because of the difficulty in separating the contribu- 
tion from highly biased intrinsically faint sources from less biased 
intrinsically brighter objects. Furthermore, existing analyses have 
shown that the simplest dusty galaxy clustering model, based on a 
linear biasi ng ansatz, is inadequate to explain small-scale cluster- 
ing power jPlanck Collaboration et al.ll201 ll ; lAddison et alj|2012j . 
hereafter A12). 

In this work, we perform a joint fit to deep number counts 
from S PIRE and the Mu ltiband Imaging Photometer for Spitzer 
(MIPS ; iRieke et aljllooi) . and well as angular power spectra cov- 
ering degree to arcminute scales from 250 to 1.4 mm (1200 to 
217 GHz) from SPIRE, Planck's High Frequency Instrum ent (HFI; 
iLamarre et al]|2010l ; |pianck HFI Core Team et alj|201 lal). the Ata- 
cama Cosmology Telescop e (ACT; Fowler et alj|2007 ; Swetz et al] 
201 ll; iDunner et alj 120121) . and the South Pole Telescope (SPT; 



Carlstrom et alj|201 ll) . We simultaneously constrain a bolometric 



dust luminosity function, the dust SED, the CIB source clustering 
properties, and the evolution of the these quantities. We require that 



our model successfully reproduces not only the clustered power 
in the power spectrum, but also the Poisson, shot-noise, power. 
Thanks to this joint analysis, we are able to demonstrate the ability 
of the angular power spectrum to provide constraints on the SEDs 
of the numerous, faint, unresolved CIB sources, in addition to de- 
veloping a model with considerable predictive power. 

Many recent dusty galaxy anisotropy analyses have inter- 
preted angular power spectra using the halo model in conjunc- 
tion with_a parametric halo occupation distribution (H.O.D.; 



I Viero et alj200g|; 


Amblard et alj201 ll; Planck Collaboration et al 


201 lllShang etal 


I2OI2I; IXia et alj|2012l; iDe Bemardis & Coorav 


2012l;IViero et alj|2012bl). In this approach, all CIB sources inhabit 



collapsed dark matter haloes, and the clustered source contribution 
to the angular power spectrum arises from correlations between ob- 
jects in the same halo ('one-halo' term, dominant on small angular 
scales) or different haloes ('two-halo' term, dominant on large an- 
gular scales). The abundance and bias of these haloes, typically 
inferred from large A'^-body simulations, is used to calculate the 
contribution of the clustered dusty sources to the angular power 
spectrum, rather than directly dealing with the dusty source bias. 
The analyses listed above make a number of assumptions, includ- 
ing; 

(i) the halo bias (relative to the linear theory dark matter power 
spectrum) may be taken as independent of scale when calculating 
the two-halo term, and 

(ii) the spatial distribution of the dusty sources within their host 
halo follows the halo mass density profile, truncated at the virial 
radius. 



ICoorav & Shethi ( I2OO2I) suggested using the scale-independent, 
large-scale halo bias with respect to the linear matter power spec- 
trum in the two-halo term as a crude way to prevent overestimating 
the two-halo power on scales where the non-linear matter power 
spectrum would begin to receive contributions from within a sin- 
gle halo (since galaxies in a single halo are supposed to be ac- 
counted for in the one-halo term). Several more recent analyses 
have found that, in fact, the halo bias with respect to the linear mat- 
ter power increases with scale for k > O.lh Mpc~ \ both at low and 
high redshift (e.g. ICole et alj|200a; Tinker et alj|2005l; 1^ et al " 



20091; iFemandez et al.ll2010l;lTinker et alj2012l;lMandelbaum et al 



2012h . Wavenumber k ~ O.lh Mpc~^ corresponds to multipole 
moment £ of several hundred at redshift 2 ~ 1 — 2, where much 
of the CIB anisotropy signal in the far-infrared and longer wave- 
lengths originates. The two-halo term dominates the one-halo on 
these angular scales, only becoming subdominant ai£ > 750, even 
when only the linear matter power spectrum is used in the two- 
halo term ^Amblard et al]|201 ll; IPlanck Collaboration et ai]l20I ll ; 



IShang et alj|2012l ; lxia et alj|2012l) . This suggests that, if the scale 
dependence of the halo bias is neglected, either power over a range 
of scales may be wrongly attributed to the one-halo term, or the 
inferred redshift distribution of the dusty sources will be altered, in 
either case likely biasing the inferred H.O.D. parameters and char- 
acteristic halo mass scales. 

The angular scale dependence of the one-halo term is deter- 
mined by the Fourier transform of the spatial distribution of sources 
within haloes (sometimes called Ugai). Setting iigai = ""dm, the 
Fourier space density profiles of a spherical NEW jNavarro et alj 
Il996l) halo, truncated at the virial radius, is a convenient choice, 
but there is no evidence that it is an accurate reflection of the 
distribution of dusty sources. There is, indeed, considerable evi- 
dence that galaxies near the center of groups and clusters, where 
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the NFW dark matter density profile peaks, are for ming stars 
less actively than tho se clo ser to the outs kirts (e.g. iKennicutll 
1983.lHashimotoetalJ 1998, iBai et al.l 2006, iBai et all 2007; also 
iBoselli & Gavazzil 2006 and references therein). These systems 
are relevant for the calculation of the one-halo power, which re- 
ceives no contribution from smaller haloes hosting a single galaxy, 
and allowing the CIB sources to lie further from the center of 
their halo than a virial radius has a non-negli gible effect on the 
one-halo power at sc ales of a few arcminutes dViero et alj|2009l : 
lAddison et ai]|2012bh . 

It should be noted that the above issues are related to how 
the halo model is implemented rather than being limitations of the 
formalism itself. It is not clear, however, that quantifying the uncer- 
tainties induced in the halo model calculations by scale-dependent 
halo bias or alternative spatial distributions of dusty sources in 
haloes may be straightforwardly accomplished. It is also unclear 
that quantities such as the halo mass function, large-scale halo 
bias, and halo concentration, are sufficiently well known at z > 2 
that their uncertainties may be safely neglected, given the excellent 
signal-to-noise of current and future SPIRE and Planck measure- 
ments. For these reasons, we do not attempt a halo model analysis 
in this work. While we are not able to directly extract (for instance) 
characteristic mass scales for haloes hosting dusty sources, we ar- 
gue that, for the reasons discussed, existing halo model analyses of 
angular power spectra in the FIR and longer wavelengths do not 
necessarily provide robust constraints on these quantities either. 

This work could be considered an extended version of the 
analysis of A12, constraining a more physically motivated model 
using higher quality spectra, as well as number c ounts. Compared 
to the CIB evolution model of iBethermin et all (2011, hereafter 
Bll), this work uses the angular power spectra and counts to- 
gether, rather than just one-point statistics (counts and LFs). Bll 
utilise data from a wider range of frequencies. Our selection is re- 
stricted to the MIPS 70 /im data and lower frequencies; while our 
SED parametrization allows phenomenologically for the presence 
of hot dust components (see Section 2.2), our focus is on the ther- 
mal dust emission associated with star formation. We also do not 
consider a model for polycyclic aromatic hydrocarbon (PAH) or 
other spectral line emis sion, which begin to contribute to the CIB a t 
higher frequencies (e.g. lLeger & Pugelll984l ; lLagache et al.l2004l) . 
An expanded SED model, and data from, for e xample, IRAS, or th e 
Wide-Field Infrared Survey Explorer (WISE; IWright et al.ll201oh . 
will be included in our model at a future date. 

The layout of this paper is as follows: in Section 2 we describe 
the model parameters we are attempting to constrain with our fit- 
ting, in Section 3 we describe our treatment of the data and un- 
certainties, results are presented in Section 4, and a discussion and 
conclusions follow in Sections 5 and 6. Calculations are performed 
using a flat, ACDM cosmology, with = 0.2715, 9.^ = 0.0 455, 
h = 0.704, Us = 0.967, and erg = 0.81 jKomatsu et alj201 ih . 



2 MODEL 

In this section we present our model for the spectral and cluster- 
ing properties of the dusty sources and show how it may be used 
to predict various statistics. Our model is divided into three parts 
- the luminosity function, which describes the abundance of dusty 
sources as a function of redshift and integrated dust luminosity, the 
spectral energy distribution, which contains the frequency depen- 
dence of the dust emission, and the clustering properties (i.e. bias). 



which connects the spatial distribution of the sources to that of the 
underlying dark matter. 

Existing studies have largely dealt with only one or two of 
these three components at a time. Bll constrained the luminos- 
ity fu nction evolution assuming fi xed source S EDs fr om earlier 
work jLagache et al.l 1 20031 12004|) . IPenin et aL I l l2012j) then cal- 
culated predictions for the angula r power spectra using the con- 
straints obtained by B 1 1 . Similarlv. lXia et aLrj2012h fixed the spec- 
tral properti es of the CIB sour ces using the results of lGranato et alj 
I2OO4}) and iLapi et alj ( [201 l|), and constrained source clustering 
properties using the spectra. We favour a combined analysis be- 
cause, while it is more complex, it is capable of obtaining more 
robust parameter constraints and performing a more thorough ex- 
plo ration of the parame ter degeneracies. Both lShangeTai] l20I2h 
and lViero et alj ( 1201 2bh attempted to constrain elements of the LF, 
SED and clustering properties simultaneously, however these anal- 
yses failed to successfully fit the Planck and SPIRE data, respec- 
tively. This may be partly because of issues with their implemen- 
tations of the halo model, discussed earlier, and also because of 
overly simplistic parametrizations of the LF and SED, and insuffi- 
cient freedom in the evolution of these quantities (discussed further 
in Section 4). 

Given that the parametrization we adopt is largely phe- 
nomenological, an important question is whether our model is ca- 
pable of reproducing any data outside that directly used to constrain 
it, and we make comparisons with several other data sets in Section 
5. 



2.1 Luminosity function 

For brevity we refer to the bolometric thermal dust luminosity, 
i'dust, integrated over 8 < A//im < 1000, simply as L. The lu- 
minosity function, "I>, is defined such that the comoving number 
density of sources with luminosity in the interval [L,L -\- dL\ is 
given by 



$(L)dlogL 



dN 



(1) 



Here, and throughout, 'log' refers to the base- 10 logarithm. We 
adopt a double-exponential form for $: 



L 



exp 



log 



1 + 



(2) 



with normalisation $c, characteristic luminosity Lc, power law in- 
dex qlp an d spread parameter ctlf. This parametrization was in- 
troduced bv lSaunders et al. ( 199 ^ ). and has been used in a range of 
CIB source analys es (e.g. IPozzietaLlbOO^; He Floc'h et alj|2005l ; 
ICaputi et al1l2007l Bll). The number of faint sources diverges for 
Qlf > 1 but we choose to allow this behaviour since, provided the 
total luminosity is convergent (i.e. qlf < 2), it is not necessarily 
prohibited by the data. 

We parametrize the evolution of the characteristic dust lumi- 
nosity, Lc, and normalisation, <l?c, using expansions about a pure 
power law in 1 + 2, 

In Lc{z) = In io + fL ln(l + 2) + Cl ln^(l + z) 
Lc{z) = Lo(l + zY^ exp [Cl ln'(l + z)] , 



(3) 



and, similarly, 

$,(2) = $0(1 + 2)'* exp [C* ln{l + 2)] . 



(4) 
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These expression reduce to a simple power law if e is the only non- 
zero evolution parameter, however more complex evolution, featur- 
ing a turning point or plateau, is allowed if the higher order parame- 
ters are non-zero. We prefer this expansion over the series of broken 
power laws adopted by B 1 1 because it ensures that all derivatives 
are continuous, and, if a quantity peaks at a particular redshift, this 
redshift will emerge naturally (provided enough terms of the series 
can be constrained), without having to be separately fitted for or put 
in by hand. The baseline model contains only as many higher order 
evolution parameters for each quantity as improve the model likeli- 
hood; some alternative or extended parametrizations are considered 
in Section 4. 

We do not find significant evidence for evolution of the shape 
parameters qlf and ctlf with redshift (see discussion in Section 
4.2). 



2.2 Spectral energy distribution 

The flux observed at frequency v from a source at redshift z with 
bolometric luminosity L is given by 



Su{L, z) 



47r(l + 2)x2(2)' 



(5) 



where x is the comoving distance to redshift z, and the differential 
luminosity, L^-, is related to the bolometric dust luminosity L and 
the source spectral function 0^ via 



(6) 



where the integral in the denominator of the fraction is taken over 
8 ^m < c/u < 1000 ^m. 

Our model is intended to phenomenologically capture the 
global properties of dust emission, such that all the parameters of 
the model are constrained from the data, while still providing a 
good fit. We consequently make several simplifying assumptions 
when modelling the dust SED. Most importantly, we consider only 
a single SED at each redshift, motivated by recent studies that have 
found relatively little scatte r about such a 'universal' SED sh ape 
for observed dusty galaxies l lElbaz et alj201 ll : lLapi et ani201 ih . 

A graybody function of the form B^{T) has been used ex- 
tensively to model FIR SEDs, where the emissivity spectral index 
/3 ^ 1 — 2 is rela ted to the emission properties of the dust grains 
( lHildebranj[l983h . Various modifications to a single-temperature 
model hav e been suggested to better fit observe d dusty galaxy 
SEDs (e.g. iBlain et alj|2003l : iHavward et alj|2012l and references 
therein) and we find that, indeed, the data reject a single-T, single- 
P dust component. We adopt a six-parameter model to describe the 
dust SED and its evolution. There are two dust components, one at 
temperature Ta, which evolves with redshift as 



Ta{z) = 



To 
To 



l+zr 



for 2 ^ 2t 
for 2 > 2t, 



(7) 



that is, constant below z — zt and then rising as a power law in 
1+2, and the other at a fixed temperature, T(,. We would expect dust 
temperature to increase naturally with incr easing redshi ft due to 
the increasing temperature of the CMB (e.g. lBlaiijll999al) . but find 
that, for the range of SED parameter values preferred by the data, 
this effect is negligible over the relevant redshift range. A single 
spectral emissivity index, /3, is common to both dust components, 
and the final parameter, /t, describes the contribution of the '6' 
dust component to the total dust luminosity. The spectral function. 



summing the contributions from both dust components, is given by 



j-du'u'PB,.{Ta) + ^'''jdu'u'^B,,{Tt,) 



(8) 



so that fb = i corresponds to the 'a' and '6' components making 
equal contributions to the total luminosity. 

We considered a range of modifications to this model, such 
as allowing a separate emissivity index for each dust component, 
or evolution in the emissivity index. For the data considered, none 
of these modifications lead to sufficient improvements in the model 
likelihood to warrant their inclusion in the baseline model (see Sec- 
tion 4.3). 

iHall et all(l201^ , IShang et"ai] ( l2012t) and lViero et all ( 1201 2bl) 



follow iBlainI ( 1 1999b ) and replace the graybody Wien tail in the 
rest-frame mid-infrared (MIR) with a power law in frequency, 
(f>iy oc as a phenomenological way to account for the pres- 
ence of hotter dust without explicitly parametrizing additional dust 
components. We choose to parametrize a second dust component 
directly in the baseline model, in order to facilitate, for instance, as- 
sessing the effect of allowing different emissivity indices between 
the two components, or evolution in their relative contribution to 
the total dust luminosity, however, we also investigated a range of 
SED models involving a power-law modification to the graybody 
Wien tail. We find that the two approaches yield similar results, 
provided evolution in dust temperature is allowed in each case (see 
Section 4.3). While an SED falling as a power law may provide 
a more realistic description of the emiss ion from very small dust 
grains shortwards of the SED peak (e.g. iDesert e t al. 1990), there 
is sufficient flexibility in our baseline parametrization to mimic this 
behaviour for the data considered. 

It is important to note that the extent to which physical dust 
properties are captured by our model is not clear Varying degrees 
of correlation between luminosity and dust temperature, not in- 



of dusty sources (e.g. iDunne et alj |2000| : 



Blain et aU 
Saiina et al] 



2003 



and references therein; 



Dunne & Eales 


2001 


Chaoman et al. 


2005 



20061 ; iGreve 67111120121) . There is, however, consider- 



able scatter in this relation, both within and between different sam- 
ples, and the extent to which it applies to the distant faint sources 
that contribute significantly to the angular power spectra, is un- 
known. Similarly, we have not considered the effect of dust self- 
absorption (the graybody equation we adopt is only valid in the 
optically thin limit), which may be significant, particularly at high 
redshift (e.g. iBenford et al.lll999l ; iBlain et al.ll2003l ; lDrainell2006l; 
iHavward et al.ll20I2h . The fact that we see good consistency be- 
tween our model predictions and various data over a range of fre- 
quency suggests our treatment is, however, at least phenomenolog- 
ically reasonable, and it seems likely that the effects of the L — T 
correlation or dust opacity, and their evolution, may be partially 
absorbed into our SED parameters (see Section 4.3 for further dis- 
cussion). 

We do not explicitly account for a duty cycle of the dust emis- 
sion in our mod el (as considered in recent fits to CIB angular powe r 
spectra by, e.g. IShang et al.ll20I2l ; lOe Bemardis & Cooravll2012h . 
We would expect our model predictions to be sensitive to whether 
star formation and dust emission on average happen continuously 
or in short-duration bursts - the Poisson anisotropy power, for in- 
stance, will be higher if, for some fixed total dust luminosity, the 
emission occurs from fewer, more luminous objects. In our model, 
this behaviour is absorbed into the luminosity function (and any 
evolution thereof), and we find excellent consistency in parameter 
constraints when fitting to the number counts alone, counts plus 
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clustered power, and counts plus clustered plus Poisson power (see 
Section 5.5), suggesting that, again, our treatment is phenomeno- 
logically reasonable. 



2.3 CIB source clustering 

We write the three-dimensional power spectrum of dusty galaxies 
as 



Pgal(fc, = {&gal(fc, Z)) PoM{k,z), 



(9) 



where the bias, (6gai), is averaged over all sources at a given red- 
shift. We model the scale and redshift dependence of the bias as 



(bgai(A:,^)> = &o(l + ^)'"' 



1 + ^bia 



exp 



^2 



(10) 



where the pivot scale, fco> is chosen to be 1 Mpc~^, and 6o, ebias, 
^bias, and kc are free parameters. In our chosen model, kc provides 
a cut-off scale at small scales, similar to that arising from tigai in 
the halo model. We choose an expansion in powers of k to mimic 
scale-dependent halo biasing, finding that the current data do not 
require a quadratic term in fc, and that its inclusion does not impact 
our results. 

A range of galaxy biasing treatments exist i n the literature; w e 
repeated our analysis using both the 'Q-model ' jCole et alj|2005h . 
with Pgai = -Pdm, and a pure power law, Pgai oc k ' 

(such that on small scales the clustered angular power spectrum 
Cf oc I '), motivated by the remarkable success of this sim- 
ple f orm in describing th e clustering of various galaxy populations 
(e.g. IWatson et al.ll201ll . and references therein). We find that nei- 
ther of these forms provides a better fit than the baseline model. 
The power law form is disfavoured at a significance level of over 
3(7, indicating, for the first time, a preference for a deviation from 
power-law clustering of unresolved CIB sources, which arises from 
the wide range of angular scale probed by SPIRE and Planck. We 
find 7 = 1.33 ± 0.03, somewhat higher than, though not inconsis- 
tent with, the value of 1.25±0.06 from A12, obtained using spectra 
spanning a similar range of frequency. 

Importantly, adopting these alternative Pgai descriptions leads 
to minimal changes in the LF and SED parameters; these other parts 
of our model are largely decoupled from the source biasing treat- 
ment due to including both counts and spectra in our fit (see Section 
4.1). 

For simplicity, we take the dark matter power spectrum in 
equation (10) as the linear dark matter power spectrum, but find 
that treating the bias as being relative to the nonlinear matter power 
spectrum instead has, again, very little effect on the LF and SED 
parameters. We calculat e the linear matter power spectrum using 
the CAMeQ distribution jLewis et al.ll2000l) . 

Having described the dusty source properties we aim to con- 
strain, we now turn to how model predictions necessary to perform 
a fit to real data may be calculated. A summary of the model param- 
eters is given, along with the marginalised parameter constraints, in 
Section 4. 



^ http://camb.info/ 
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2.4 Angular power spectrum 

The angular power spectrum of CIB anisotropics from correlating 
a pair of maps at frequency v is written as the sum of cl ustered and 
Poissonian components (e.g. IPeeblesll 1 98d ; ISondl 1 996h : 



(11) 



where I is the multipole moment. We take the Poisson contribu- 
tion, which can be thought of arising from the correlation of the 
image of a source in one map with the image of the same source in 
the other map, to be independent of angular scale. On scales larger 
than 18" (the beam FWHM of the 250 /im SPIRE detector), all 
but fairly nearby galaxies will appear as point sources, and we as- 
sume that any local, extended sources, which may contribute scale- 
dependent Poisson power, since some internal structure is resolved, 
are masked from maps prior to calculating spectra (for instance by 

cross-matching with existing catalogu es). 

In the flat-sky limit jLimbej Il953l ; iKaiseJ Il992h . the 
clus tered power is written as a single line-of-sight integral 



(e.g.lBond et al.lll991al;lKashlinskv & Odenwaldll2000l ; lKnox et al] 



I2OOI ; Tegmark et alj 2004). Since the number of faint sources is 
not required to be finite in our model, it proves helpful to work 
in terms of the comovi ng emissivity density, jllaiman &Knoxl 
l200(]| ; lKnox et al.l200ll) . If source bias and luminosity are uncorre- 
lated (see below), the clustered power can be written as 



dz dx 1 



{b^^x{k,z)) (>(2))^mPDM(fc,z), 



(12) 



where wavenumber fc = (£ + 1/2) /x[z), and (ji.)cut is the mean 
emissivity density once bright sources (with flux 5* > S'cut) have 
been masked from the map, given, in terms of the quantities defined 
above, by 



(>(2)>cut = (1 + Z)X 



dL HL,z) 
"dS,, Lin 10 



(13) 



where here L is considered a function of Sv, by inverting equations 
(6) and (7). The fact that Sv and L are related by a simple constant 
at ea ch redshift in our model simplifies this calculation. 

iBethermin et al.l l l2012cl) recently found evidence for depen- 
dence of star formation rate (and thus bolome tric IR luminosity) 
on halo mass through abundance matching, and lWang et alj ( 1201 2h 
also found that the SFR-halo mass relation i s not flat but has a peak 
at some characteristic mass scale (see also lSehroozi et aLll2012bl 
and references therein). Since the halo bias i ncreases with mass 
(e.g. lSheth & Tormeij[l999l ; [Tinker et al.ll2010t) . we may therefore 
expect a correlation between dust luminosity and bias, and thus be- 
tween j,j and fogai. To allow for this, we consider re- writing the 
clustered power, via 



(&gal(fc,z))(>(2:))cut -> (fegal(fc,2)>(z))cut. 



(14) 



We investigate a range of treatments for the mean product of fcgai 
and ji,, and find that, for the data considered, allowing a 6gai — L 
correlation serves primarily to degrade constraints on 60, and does 
not have a significant effect on either the goodness-of-fit or the LF 
and SED parameters. We therefore do not include such a correlation 
in the baseline model; more discussion is provided in Section 4.4. 

2.5 Poisson anisotropy power 

Existing studies have typically treated the Poisson power either 
by adding an additional free parameter to each spectrum (e.g. 
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lHaiianetal.ll2012l : lAmblard et al] |201 1': 'Xia et al."2012'), or using 
the predictions of the B 1 1 CIB model ( Planck Collabor ation et alj 
l201lUShangetai]|2012h . We instead require that our model for the 
dusty source LF and SEDs successfully reproduce the Poisson, as 
well as the clustered component of the angular power spectrum. In 
terms of the quantities introduced earlier, the Poisson contribution 
to the anisotropy power is given by 



, dx 2 



r 

Jo 



' dL ^{L,z) ^2 
'^"dS. Lin 10 



(15) 



Since the Poisson power is weighted more towards the brighter 
sources, one may expect that fitting to deep number counts already 
constrains the Poisson power We find that, in fact, this is not the 
case, and that the Poisson tail in the SPIRE, ACT, and SPT spec- 
tra carries additional constraining power (see Section 5.5). The fact 
that there is useful signal in the shot-noise is one of the main results 
of this work. 

Removing the contribution from sources above Scut in the in- 
tegrals in equations (13) and (15) does not exactly mimic the re- 
moval of bright sources in confused maps, which typically involves 
masking a region of roughly a beam size per source, and may there- 
fore also remove flux from from objects in the same group or clus- 
ter In Section 5.4 we show that any bias introduced by this issue 
is likely to be small. Constructing a simulated CIB sky map from 
our theoretical model, and utilising the actual masking scheme that 
was used on the real data, will be used as a more rigorous test in 
future work. 



2.6 Differential number counts 

The differential number counts are modelled using the quantities 
introduced above as (see, e.g. discussion in Bl 1) 



dN 



dx 2 dL $(L,2) 
^ dz^ dS^ LlnlO ' 



(16) 



where, as in equation (14), L is taken as a function of 5*1, and z. 
The Euclidean-nomalised counts (multiplied by a factor of S^'^) 
are integrated over the finite width of each flux bin to compare with 
the measured counts. 



3 DATA TREATMENT AND MODEL FITTING 
3.1 Data sets used 

We constrain the model parameters described in Section 2 using a 
joint fit to the following data sets: 

(i) differential number cou nts from 5pi?zer-M IPS at 70 and 160 
^m (4300 and 1900 GHz; Be thermin et"ai]|2010l) 

(ii) differential number counts from / /erjc/;e/-SPIRE at 250 , 
350 and 500 ^m (1200, 857 and 600 GHz jBethermin et alj|2012d. 
hereafter B 1 2) 

(iii) angular power spectra bandpowers from SPIRE at 250, 350 
and 500 ^m (1200, 857 and 600 GHz; Amblard et al. 201 1) 

(iv) angular power s pectra bandpowers from Planck-W Fl at 857, 
545, 353 and 217 GHz jPlanck Collaboration et al1l20I ll) 

(v) angular p ower spectrum bandpowers from ACT at 218 GHz 
( lDasetal.ll20Ilh 

(vi) angular power spectrum bandpowers from SPT at 220 GHz 



( iReichardt 



ular power i 
et al]|20I2h 



The MIPS numbe r counts are described in Tables 3 to 6 of 
ISethermin et al.l j20I(t) . and the SPIRE number counts in Tables 2, 
3 and 4 of BI2. The baseline model does not include the SPIRE 
counts obtained by stacking from the GOODS field (see Section 
3.9); we also exclude the highest flux bin in the MIPS 160 ^.m 
stacked counts since the lowest flux bin from the MIPS 160 /im 
resolved counts covers the same flux range with a smaller uncer- 
tainty. 

We summarise key properties of the power spectra in Table 2. 
Note that we do not use the full range of SPIRE or ACT bandpow- 
ers in the baseline fit; we use only data from I > 2000 for SPIRE, 
and I > 2400 for ACT, because of uncertainty over Galactic cirrus 
contamination and spectrum-to-spectrum conelations in the SPIRE 
data, and because data from lower £ in ACT do not contribute sig- 
nificant constraining power once the primary CMB power spectrum 
is subtracted (see Sections 3.8.2 and 3.9, below). Bright sources in 
the ACT and SPT maps were masked based on flux at the 150 GHz 
bands of these instruments, not 220 GHz. Our model allows us to 
account for this, since, for given model parameters, we can predict 
both the 150 and 220 GHz flux of all sources, and so calculate the 
equivalent Scut at 220 GHz and substitute these values into equa- 
tions (13) and (15). 

We do not consider counts from lower wavelengths than the 70 
fim MIPS band because, as stated earlier, our focus is on thermal 
dust emission rather than other contributions - such as emission as- 
sociated with AGN, whic h may be non-negligibl e in, for instance, 
Spitzer 24 fim counts (e.g. lBethermin et alj20I2an . or PAH features 
- to the SED. Counts from the Hersche l Photodetector Array Cam- 
era and Spectrometer teertaet alj20Ilh are not included as they are 
in good agreement with the MIPS counts and span a smaller range 
in flux. We do not use counts from higher frequencies than 600 GHz 
because of possible enhancement of the intrinsic counts by strong 
gravitational tensing (see Section 3.10). Our choice of spectra was 
limited lou ^ 220 G Hz due to the presence of the th ermal Sunyaev 
Zel'dovich effect (tSZ; ISunvaev & Zerdovicljl970l) , possible tSZ- 
CIB c ross-correlations ( e.g. IShirokoff et al]|201 1 ; IReichardt et alj 
I2OI2I ; IZahnet all [2012! ; lAddison et al.ll20I2bl) . and the increased 
importance of the primary CMB anisotropics, in data at lower fre- 
quencies. 

Towards the end of this study, IViero et si\ l l20I2bl ) presented 
new SPIRE power spectrum results. We will consider these data 
in future work; for the p urposes of this paper we n ote that the 
SPIRE bandpowers from lAmblard et alj 1 I2OI ll) and IViero et al] 
( l2012bh are, for £ > 2000, in agreement within a multiplicative 
shift of 10 per cent (Figure 7 of lViero et aLll20I2bl) . which is small 
compared to the flux calibration uncertainty of the maps used by 
lAmblard etal] ( I20T l|) of 15 per cent (corresponding to a > 30 per 
cent uncertainty in units of power). 

3.2 Model likelihood and data covariance 

We explore the par ameter space using a M arkov Chain Monte Carlo 
analysis (MCMC; iMetropolis et aljl 19531) using the optimal sam- 
pling step size from lDunklev et aLr ( l2005l) - see also lGelman et alj 
( Il996l) . The results presented later in this paper were obtained from 
chains of approximately 6 x 10^ steps. Running longer chains 
did not lead to significant shifts in the parameter covariance or 
marginalised constraints. Evaluating the full posterior probability 
at each step can be achieved in under a millisecond on a modern 16- 
core processor, meaning even these long chains take no more than a 
day to run. The long chain length required to achieve convergence 
is a consequence of the complex parameter degeneracies rather than 
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Table 1. Angular power spectra used for constraining tlie baseline model 



Instrument 


1^0 

(GHz) 


Bandpowers 


Angular' scale 


91" 

•^cut 

(mJy) 


(per cent) 


Reference 


Herschel-SVmE 


1200 


23 


2000 < 1 


? < 63300 


50 


15 


Amblardetal. (2011) 




857 


21 


2000 < ( 


? < 41200 


50 


15 




600 


19 


2100 < t 


? < 26300 


50 


15 




Planck-HFI 


857 


9 


80 <e 


< 2240 


710 


9* 


Planck Collaboration et al. (2011) 




545 


9 


80 <e 


< 2240 


540 


9* 




353 


9 


80 <i 


< 2240 


325 


4* 






217 


9 


80 <e 


< 2240 


245 


4* 




ACT 


219.6 


13 


2391 < 


e < 9900 


20t 


7 


Das etal. (2011) 


SPT 


219.6 


15 


2000 < 


e < 9400 


6.4t 


4.8* 


Reichardt et al. (2012) 



* level above which bright sources are masked in the map prior to calculating spectra - Scut is given at 150 GHz for ACT and SPT (see text) 
* absolute map calibration uncertainty: given in units of power for the SPT spectrum, otherwise in units of flux or temperature 
we take the Planck calibration uncertainties to be 9 or 4 per cent, rather than the nominal 7 or 2 per cent, to account for uncertainty in the bandpass filter 

responses (Section 3.6.1) 



solely the number of parameters (which - including model param- 
eters and the additional parameters described below - is 40 for the 
baseline model). 

The posterior probability is proportional to the product of the 
likelihood, £, and the prior probability. We work with the negative 
log-likelihood, given by 



-21n£ = Y^[M,{9) 



[M>(6() -D,] 



(17) 



where i labels the different data sets, the elements of the data vec- 
tors Di are either the binned spectrum bandpowers (see below) or 
the binned number counts, as appropriate, Mi(S) are the corre- 
sponding vectors of model predictions, which depend on the model 
parameter values at each step, denoted by 6, and Ct are the covari- 
ance matrices, which contains the data uncertainties. 

Note that we include off-diagonal data covariance elements 
for many of the data sets. These are discussed further in Sections 
3.7 and 3.9. 

At each step of the MCMC chain, the likelihood is multiplied 
by the contribution from parameter priors described later in this 
section. In particular, Gaussian priors are adopted on the photomet- 
ric calibration parameters, /cai, which are discussed in Section 3.5. 

The remainder of this section contains discussion of vari- 
ous issues relating to fitting models to the counts and angular 
power spectra. We draw the reader's attention to our treatment of 
the tension between the SPIRE a nd Planck spect r a discu ssed by 
IPlanck Collaboration et akl ( 1201 ih and I Viero et"al] ( l2012bl) . These 
papers suggest that the discrepancy may arise from incorrect cal- 
culation of the SPIRE effective beam areas, or uncertainty in the 
Planck beam shape, respectively, and we allow for both these pos- 
sibilities as described in Section 3.7, below. 



3.3 Binning of spectra in £ 

The measured power spectra were calculated using bins in £ rather 
than individual £ values. The resulting bandpower values Ct are 



related to the unbinned Cis by 



(18) 



where the sum is over all £ values in bin 6, and Wb,e is a weight 
function, which is equal to unity for th e SPIRE spectra, and I for 
the P lanck spectra (see Section 4.1 of IPlanck Collaboration et al] 
l20Ilh . The weighting of the A CT and SPT spectra is slig htly more 
complex; see lDas et al.l ( 1201 ih and lReichardt et al] ( 1201 2h for more 
details. We use the same binning and weighting scheme used for 
the measured bandpowers when calculating the model predictions. 



3.4 Correlated uncertainties for staclced counts 

The deepest MIPS and SPIRE counts used in our fit were obtained 
by stacking on 24 /im Spitzer sources. This process is subject 
to systematic uncertainties that may be expected to b e correlated 
across bins in flux and also from band to band (e.g. Iviero et al] 
l20I2ah . This is supported by the fact that the scatter in the stacked 
counts is smaller than the quoted uncertainties. We assume a 50 
per cent covariance between the uncertainties in the stacked MIPS 
counts at 70 fim and 160 pim, and a 50 per cent covariance between 
the uncertainties in the stacked SPIRE counts between all flux bins 
and across all three bands within each redshift bin. While our re- 
sults are largely unaffected by changes in the assumed covariance 
by several tens of per cent, more detailed analysis of the covariance 
in future deep number count extraction is important for ensuring 
robust constraints can be obtained. 



3.5 Photometric calibration uncertainties 

The data considered in this work were obtained using in- 
struments that do not measure absolute flux and must 
therefore be c alibrated using ancill ary meas urements (see 
IStansberrv et alj 2007. [Gordon et all 2007. ISwinvard et al] 



2010. IPlanck HFI Core Team et alj 201 1, iHaiian et al 



2011, and 



iLueker et al.l 2010. for details of the MIPS 160 pm, MIPS 70 fim. 
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SPIRE, HFI, ACT, and SPT data, respectively). Each power spec- 
trum was calculated from maps with the best-guess flux calibration 
correction applied, however the uncertainty in this calibration 
is not negligible compared to the statistical uncertainties in the 
bandpowers, and it is therefore necessary to correctly account for 
the calibration uncertainty during model fitting. 

Let the calibration of the maps for the quoted bandpower val- 
ues equal unity, and let /cai be the factor, in flux units, that the map 
must be multiplied by in order to achieve the correct calibration. We 
fit for /cai as nine extra free parameters (one for each spectrum), as 
in A12. At each step of the MCMC chain, the bandpowers and er- 
rors of each spectrum are multiplied by the corresponding /^ai be- 
fore the likelihood is calculated. Furthermore, in order to penalise 
models requiring /cai to be considerably different from one, we im- 
pose Gaussian priors on each /cai, with mean of unity and standard 
deviation equal to the quoted calibration uncertainty for that band. 
We also account for an additional effect of the calibration uncer- 
tainty on the power spectrum, namely that bright sources are not 
masked to the nominal flux value. Scut, but actually to /cai Scut, by 
substituting this modified value into equations (13) and (15). 

Absolute flux calibration uncertainty also affects the differen- 
tial number counts; the lower and upper limits of each flux bin, 
Smin and Smax, must be replaced by /caiSmin and /caiSmax, re- 
spectively, at each MCMC step, while S^''^^ - > fcafS ^ '^f^- In 
the time between the analyses of lAmblard et al.l ( l201ll ) and B12, 
the SPIRE flux calibration uncertainty was improved from the 15 
per cent quoted above to 7 per cent in each band, with a band- 
to-band relative calibration uncertainty of 2 per cent (values taken 
from the SPIRE Observers' ManuaQ, hereafter SOM). We account 
for this by fitting for separate calibration parameters for each band 
for the SPIRE spectra and the SPIRE counts. Two additional cal- 
ibration parameters are included for the MIPS 70 and 160 /im 
counts. 



3.6 Accounting for bandpass filter profiles 

The CIB intensity and anisotropy depend str ongly on frequency 
over the wavelength ranee con sidered (e.g. iFixsen et alj 119981 ; 
IPlanck Collaboration et a"lll201 ih . with the anisotropy power at a 
given angular scale increasing by up to an order of magnitude 
across the microwave Planck filters (see Figure 3 of A 12). Ac- 
counting for the bandpass filter transmission profile when fitting 
models to the measured data is therefore important. We follow 
IPlanck Collaboration et al.l | |20I ih and A12 in working with the in- 
tegrated CIB SED rather than calculating filter responses to indi- 
vidual source SEDs at all redshifts. 



3. 6. 1 Planck spectra 

We convert the Planck bandpowers from temperature to flux den- 
sity units using the values given in th e notes accompanying Table 
4 of IPlanck Collaboration et al] ( 1201 ll) . This conversion assumes a 
source SED that varies as S,^ oc 1/v (i.e. vS^ = const.). Convert- 
ing the spectra to the 'real' flux units for a source with observed 
SED requires multiplying both the bandpowers and errors by a 



^ http://herschel.esac. esa.int/Docs-SPIRE/htmiyspire_om.html#xl- 
880005.2.1 



factor /sED, where (e.g. lPlanck HFI Core Team et alfcOl Ibl) 



/sED 



(19) 



with bandpass filter transmission profile t{i') and nominal 

band frequency vp (the se values are listed in Table 1). 

IPlanck Collaboration et al. I boiih report 1//sed values of 1.00, 
1.06, 1.08, and 1.08 at 857, 545, 353 and 217 GHz, respectively, 
calculated using th e CIB SED measured from FIRAS data by 
iGispert et af] ilOOdi . We find that the CIB SED predicted by our 
best-fit model is very similar in shape to that measured by FIRAS 
(see Section 5.1), and do not make further corrections to the Planck 
bandpowers. Uncertainties in the Planck /sed val ues are estimated 
at 2 per cent jPlanck HFI Core Team et al .1201 id) ; to allow for this 
we take the Planck calibration uncertainties to be 9 per cent for the 
857 and 545 GHz bands, and 4 per cent for the 353 and 217 GHz 
bands, rather than the nominal 7 and 2 per cent. 



3.6.2 SPIRE spectra and counts 

The SPIRE maps were made assuming point source emission and 
(X l/u (see the SOM for further information). For the SPIRE 
spectra, we calculate /sed values using Table 5.3 of the SOM for 
extended source emission and input SEDs of the form S 1/ oc f 
with effective spectral indices, a, of 0.0, 1.0, and 2.0, for the 
1200, 857 and 600 GHz bands, respectively. These values were se- 
lected based on the frequency dependence of [oispert et al. ( 20()3) 
CIB SED mentioned above; the best-fit CIB SED from our model 
yields very similar values. A separate factor, denoted Ka,eIK^p 
in the SOM, is also applied to account for the detector response 
to extended as opposed to point-like emission. These factors to- 
gether yield effective /|ed correction factors of 0.989^, 0.991^, 
and 0.995^ for the 1200, 857 and 600 GHz SPIRE spectra, respec- 
tively. 

For the SPIRE counts, the SED correction is applied as 5* — >^ 
/sed5, and S^-^dNjdS fsioS^-^dN/dS, where here the 
conversion factors are for point-like emission, and we do not in- 
clude the KiE/Kip connection, yielding /sed values of 0.989, 
0.974, and 0.940 at 1200, 857 and 600 GHz, respectively. 



3.6.3 ACT and SPT spectra 

The ACT and SPT bandpowers are reported in units of CMB tem- 
perature, which can be converted to flux density (i.e. pK? to Jy^^ 
sr^^) by multiplying by (c'i3^/9r|T=TcMB where 



dT 



c2 (e- - 1)= 



(20) 



X — h^v /k^TcMB, and Tcmb is the present-day CMB temper- 
ature. We adopt an effective frequency, i^o, of 219.6 GHz for per- 
forming this units conversion and evaluating the model predictions 
for the ACT spectra; this value accounts for the ACT bandpass filter 
and assumes an SED rising as Sv oc v^'^ jSwetz et al.ll201 ih . We 
assume that small differences in SED or uncertainties in the effec- 
tive frequency can be absorbed into the ACT spectrum calibration 

unce rtainty. 

iReichardt et alj | |20I2|) report the same effective frequency as 
ACT, 219.6 GHz, for the response of the SPT 220 GHz channel to 
sources with S,^ oc v'''''\ and we use this value for calculating the 
units conversion and model predictions for the SPT bandpowers. 
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3.6.4 MIPS counts 

iBethermin et alj ^20 id) reported MIPS counts assuming, again, a 
source SED that varies as oc 1/v. We calc ulate the MIPS /sed 
values by integrating the lGispert et alJ ( 120001) CIB SED across the 
MIPS bandpass filters (equation 19), finding 0.977 and 1.012 at 70 
and 160 /im, respectively. We assume, as above, that small uncer- 
tainties in these values can be absorbed into the flux calibration 
factors. 



3.7 Instrumental beam treatment 

Given t he observed discrepancy between SPIRE and Pl anck 
spectra jPlanck Collaboration et alj|201 ll ; IViero et alj|2012bl) . we 
opt to allow additional freedom i n the SPIRE and P lanck 
beam treatment compared to the lAmblard et alj j201lh and 
IPlanc k Collaboration et aljj 201 ih anal vses. 

|pianck Collaboratio n et alj ilOl ll) found from a preliminary 
SPIRE / Planck cross-correlation that the discrepancy between the 
SPIRE and Planck spectra could be resolved if the SPIRE beam ar- 
eas were underestimated by around 5 and 10 per cent at 857 and 600 
GHz, respectively. We therefore introduce an additional beam area 
correction parameter, /beam, for each SPIRE band, with a Gaussian 
prior centred at unity with a Icr uncertainty of 10 per cent. At each 
MCMC step the SPIRE bandpowers and errors are multiplied by 



the corresponding /beam before the likelihood is calculated. 

For the Planck spectra, we use the beam error estimates given 
in Table 4 of[pianck Collaboration et al. (201 1), however we mul- 
tiply the Planck beam errors by five before taking the quadrature 
sum of the statistical and beam uncertainties, motivated by the pos- 
sible SPIRE and Pl anck spectrum shape discrepancy discussed by 
IVieroetaljj2012d) . We also assume a 100 per cent correlation in 
the beam uncertainty across bandpowers in each Planck spectrum, 
appropriate for a roughly Gaussian beam whose width is uncertain. 

Allowing a larger SPIRE beam area correction, or increas- 
ing the Planck beam shape uncertainties further, ha s no impact on 
our re sults. If the nominal treatments present ed in lAmblard et alj 
( l201lh and IPlanck Collaboration et alj ( 1201 ih are used, the mis- 
match between the SPIRE and Planck spectra is largely absorbed 
by the SPIRE calibration parameters, with shifts in the CIB model 
parameter constraints by no more than O.Sa; our choice of beam 
treatment does not have a significant impact on our results (see Sec- 
tion 4.5). 

We use the ba ndpower covariance matrix provided by 
iReichardt etal] ( |2012|) . which includes beam uncertainties, when 
calc ulating the SPT spectr um l ikelihood contrib ution. 

lAmblard et al] ( 1201 ll) and lPas et akl ( 1201 ll) found that bin-to- 
bin correlations in the SPIRE and ACT spectra, due either to corre- 
lated beam shape uncertainties, or the source mask, are small, and 
we also therefore neglect them in this work. 



3.8 Non-CIB contributions to spectra 

All the spectra used in this work were calculated from regions of 
the sky that are known to be relatively free from Galactic dust con- 
tamination, however diffuse emission (cirrus) is present; we also 
consider several additional non-CIB contributions to the microwave 
ACT and SPT spectra. 



3. 8. 1 Planck spectra 

The cirrus was cleaned from the Planck maps using cross- 
correlations with high-resolution Galactic HI maps (Section 2.5 of 
IPlanck Collaboration et alj|201ll) . and the 143 GHz Planck chan- 
nel was used to remove the primary CMB contribution. We do not 
consider any modifications to these methods in this analysis. 



3.8.2 SPIRE spectra 

lAmblard et al. l|) remove the cirrus in the SPIRE power spec- 
trum using an extrapolation fro m lower wavelengths, however 
IPlanck Collaboration et alj ( 1201 ih found that this treatment over- 
estimated th e cirrus contamination. We the refore adopt a similar 
treatment to IPe Bemardis & Cooravl ( |2012|) . and fit for the cirrus 
contamination in the SPIRE power spectra along with the CIB 
model parameters. We further reduce the possible impact of in- 
correct cirrus modelling by only including SPIRE bandpowers at 
£ > 2000 in our fit. We assume that the cirrus power spectrum has 
the form Cf" = Ae»rr(^/20 00)~""" (e.g. iGautier et alJI 19921: 
iMiville-Deschenes et ani2007l) . We fit for the amplitude Adw at 
1200 GHz, and assume a graybody frequency dependence of the 
cirrus flux with emissivity index /3 = 1.5 and e ffective temperature 
T = 20 K (motivated by the findings of e.g. 'Bracco et al j|201ll) 



to scal e the cirrus to 857 and 600 GHz. Planck Collaborati on et alj 
( I20T l|) found that the cirrus at these lower frequencies is negligible 
in the field from which the SPIRE spectra were calculated, and so 
we do not expect our results to be sensitive to the exact frequency 
dependence adopted. We adopt a Gaussian prior on ricirr, centred 
at 2.89, with a la sp read of 0.44, using the constraint obtained by 
iLagache et alj ( |2007|) . but doubling the width of the uncertainty to 
account for possible variation of the cirrus index with frequency. 
We find that, in fact, even the cirrus in the 1200 GHz SPIRE power 
spectrum is not significantly detected in our fit, and that the ex- 
act treatment of the SPIRE cirrus has negligible impact on the CIB 
constraints. 



3.8.3 ACT spectrum 

We include in our model a lensed CMB component with a shape 
in £ calculated assuming our fiducial cosmology using CAMB, but 
fit for the an overall amplitude, Acmb, as a free parameter, using 
a Gaussian prior with a Icr uncertainty of 10 per cent. This free- 
dom is more than sufficient to mimic changing, for example, Js by 
±0.025 (roughly the la constraint obtained bv lKeisleretalj201ll) 
on the angular scales on which the CMB makes a non-negligible 
contribution to the total spectrum. A more detailed treatment of the 
CMB power spectrum has no impact on our results, at least for a 
standard cosmological model. 

No subtr action of Galact ic dust is performed for the ACT 
spectrum (see iDas et alj|201ll) . We include the amplitude of the 
subdominant ACT radio Poisson component as a nuisance param- 
eter with a Gaussian prior of 4.1 ± 0.8 p.K ^ (or 0.48 ± . 09 ly ^ 
sr^ ), adopting the mean value obtained bv lDunklev etalj ( l201ll) . 
but conservatively doubling the measured uncertainty. 

While the ACT 218 GHz band is virtually at the thermal 
Sunyaev Zel'dovich effect null, the blackbody kinema tic Sun- 
yaev Zel'dovich effect (kSZ; ISunvaev & ZeldovicljlT980h will be 
present. We allow for a non-zero kSZ component, fitting for the 
power at £ = 3000 in units of £{£ + l)Ce/2n, ^ksz, using a 
fixed shape in £ obtained fr om recent hydrodynamical simulations 
dBattagha et"ani2010ll2012l) . and imposing a uniform prior of < 
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A\sszl /jK'^ < 8, based on constr aints obtained bv lDunklev et alj (assuming the same binning of the two spectra), by 



( I2OI II) and lReichardtetalJ ( l2012h . In the latter work it was found 
that a larger kSZ amplitude was possible only for very high degrees 
of correlation between the tSZ and CIB (over 50 per cent in units 
of power at i = 3000), f or which there appears little physical mo- 
tivation (see discussion in lAddison et a"l]|2012bl) . The choice of the 
kSZ amplitude prior does not have a significant effect on any of the 
CIB model parameters. 

In principle, power spectrum measurements at ~220 GHz 
could allow meaningful kSZ constraints to be obtained in the ab- 
sence of the tSZ and tSZ-CIB contributions. We find that, for the 
data considered, the dominant dusty source contribution is not suf- 
ficiently well-constrained for the kSZ constraints to be useful (see 
Section 5.3). 



3. 8. 4 SPT spectrum 

The lensed CMB and kSZ power are included when modelling the 
SPT spectrum in the same way as for ACT. Radio source Poisson 
power is minimal at 220 GHz with the SPT source mask; we in- 
clude a radio Poisson c omponent with amplitu de fixed to the mean 



of the prior adopted bv lReichardt et al 



j2012l). which corresponds 
to 0.16 Jy^ sr~^ in flux density units. iReichardt et alj j2012h also 
find that Galactic cirrus contributes only at the per cent level to the 
total 220 GHz spectrum, and we adopt the model described in their 
Section 5.6, without adding any additional free parameters. Allow- 
ing more freedom in the SPT radio source or cirrus levels has a 
negligible impact on our results. 



3.9 Cosmic and sample variance 

The deepest BI2 number counts were obtained from an area of 
only half a square degree (GOODS-N field). It is apparent from 
Figure 10 of B 12 that the counts, binned by redshift, are not al- 
ways consistent between GOODS and the larger COSMOS field 
within the quoted error bars. There is good consistency between 
the counts in the l owest GOODS flux bin and those found from the 
P{D) analysis of iGlenn et al] bOlOl) , but, since both analyses use 
the same field, this does not help inform our treatment of possible 
cosmic variance uncertainty. We choose to include only the counts 
from the larger COSMOS field counts in our fit. While we can- 
not rule out the possibility that the faint counts in the smaller field 
are, in fact, more representative of the global population, we find 
no current motivation fo r favouring these data; we also note that 
iBethermin et alj l l20I2ah recently found that the GOODS counts 
systematically fell below their CIB galaxy model predictions while 
the COSMOS counts were fairly well-reproduced. We consider the 
effect of including the GOODS counts in Section 4. 

On large scales, uncertainties in the Planck and SPIRE band- 
powers are dominated by the limited sampling of Fourier modes, 
due to finite sky coverage. For binned bandpower, Cb, with bin b 
spanning ^min ^ ^ ^ ^ max, this cosmic variance contributi on is 
given analytically by (e.g. lFowler et alj[2010l ; lDas et alj|201lh 



e ■ 

mm 



sky 



(21) 



where /sky is the fraction of sky covered by the map used to cal- 
culate the spectrum. For spectra calculated from the same patch of 
sky, there will be a correlation in this cosmic variance uncertainty 
between spectra at frequencies v\ and V2 , given, for bandpower b 



/sky 



(22) 



where j^^^, is the binned cross-spectrum. The effect of the 
spectrum-to-spectrum correlation in the SPIRE spectra may be ex- 
pected to be larger than for Planck, because around ten times less 
sky was used. Due to this fact, and the uncertain cirrus contami- 
nation on large scales in the SPIRE spectra (Section 3.8.2), we do 
not include any SPIRE bandpowers from £ < 2000, preferring the 
large-scale power to be constrained by Planck alone. Removal of 
the large-scale SPIRE data does not significantly degrade parame- 
ter constraints. 

We find that we can adequately account for the Planck 
and remaining SPIRE spectrum-to-spectrum correlation with a 
simplified iterative approach (similar to that described in A12), 
without requiring additional calculations at each MCMC step. 
We re-run the MCMC chain several times. At the end of each 
chain, we use the best-fit model parameters to calculate the ra- 
tio Cb,,,-^,,^/ {Cb,^iCb,v2) for each bandpower and for each cross- 
spectrum (e.g. 857 X 545, 857 x 353, 857 x 217, 545 x 353, 
545 x 217 and 353 x 217 GHz for Planck, where the multiplication 
symbol here denotes cross-correlation). This ratio is then multiplied 
by the corresponding measured auto-spectra bandpowers to calcu- 
late an estimate for the cross-spectrum bandpowers, which are used 
to calculate the off-diagonal covariance elements, crf ^^^, , for use 
in the subsequent chain. We find that good convergence is achieved 
after three chains, and that there is good agreement between the 
SPIRE cross -correlation ratios o btained by this method and those 
measured bv lVieroetai] ( l2012bl) . discussed further in Section 5.5. 

We note that there is a partial overlap be tween the Lockman 
Hole region used for the lAmblard et al.l ( l20Ilh SPIRE analysis and 
one of the regions used to calculate the Planck spectra. Since five 
separate, larger, regions are also used in the Planck analysis, we 
neglect the effect of any correlation in the cosmic variance uncer- 
tainty between the Planck and SPIRE spectra. There is also partial 
overlap between the ACT and SPT coverage, but as we are fitting 
to these spectra on angular scales where noise dominates the error 
budget, we similarly neglect any correlation between the ACT and 
SPT spectra. 



3.10 Effect of strong gravitational lensing 

The observed bright-end CIB source counts in the submm and 
at longer wavelengths are expected to differ from the intrinsic 
counts due to the effect of stro ng gravitational lensing by fore- 
ground groups and c l usters (e.g. |Blainlll996l; IPerrotta et al I2OO2I ; 
iNegrello et al] |2007| ; iLimaet alJ l2010bl ; iNegrello et al.l UOld) . 
While, for any given source, this effect is not frequency dependent, 
the fact that lower frequencies receive a larger relative contribution 
from sources at higher redshift, where the lensing cross-section is 
higher, means that we may expect the enhancement of the intrinsic 
counts to increase with decreasing frequency. 

A full treatment for the effect of lensing requires mod- 
ellin g the distribution o f lens systems and lensing cross-sections 
(e.g. lUma et alj l2010d ; iHezaveh & HoldeiJ 1201 ih . Furthermore, 
iMead et all | |2010|)~ find that baryonic physics, including active 
galactic nuclei feedback, may d ouble the strong lens ing cross- 
section for massive clusters, while iHezaveh et ai]i l20I2h show that 
compact sources are preferentially more strongly lensed, introduc- 
ing a source-size dependence. Accounting for the lensing effect on 
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CIB counts thus requires extensive modelling beyond that used to 
describe the intrinsic AN /AS. 

In this work, we conservatively choose not to constrain our 
model with counts from shorter frequencies than the SPIRE 600 
GHz band. Even at this frequency, signi ficant enhancement of the 
counts is p ossible for 5* > 100 mJv (e g.|Negrello et al.l2007l B 1 1, 
IWardlow et al. 2012i . lBethermin et alj|2012ah . The uncertainties in 
the B12 counts at the bright end are large due to the limited sky 
coverage, suggesting that the bias introduced in our parameters by 
ignoring the effects of the lensing is likely to be small. We consider 
the effect of excluding either the brightest SPIRE counts (5* > 60 
mjy), or the high-redshift (z 2) SPIRE counts completely, and 
find that, indeed, there is no systematic change in the parameter 
constraints. We ignore the effect of strongly-lensed sources on the 
angular power spectra because the lensing preserves surface bright- 
ness, and it is fluctuations in surface brightness that the spectrum 
measures. In addition, the brightest sources are masked from the 
maps before the spectra are calculated. 



4 RESULTS 

4.1 Goodness of fit and parameter constraints 

Globally, the model provides a good fit to the data, with 
X^/d.o.f. — 239/248 = 0.96. The degrees of freedom are cal- 
culated by summing the number of bandpowers (36 Planck, 63 
SPIRE, 13 ACT and 15 SPT) and counts (105 SPIRE and 39 Spitzer 
flux bins), and subtracting the number of free parameters (18 CIB 
model parameters, plus the SPIRE cirrus amplitude and index, and 
kSZ, primary CMB and ACT radio source amplitudes). 

The SPIRE, Planck, ACT and SPT spectra are individually 
well-fit by the model, with xVd-o.f. of 54/63, 35/36, 11/13, and 
14/15, respectively. Likewise for the SPIRE and MIPS counts, with 
X^/d.o.f. of 80/105 and 34/39. If the analysis is performed with- 
out allowing the SPIRE beam area correction or increased Planck 
beam shape uncertainty discussed in Section 3.7 the model remains 
a good fit, with x^/d.o.f. = 1.04. 

The one-dimensional marginalised parameter constraints, and 
CIB parameter correlation matrix, are shown in Tables 2 and 3, re- 
spectively. The parameters are grouped into three sections - lumi- 
nosity function, SED and bias. While there are strong correlations 
within each section, a combined fit to the counts and spectra greatly 
reduces the correlation between parameters in different sections. In 
particular, the parameters describing the CIB source bias parame- 
ters are virtually decoupled from the LF and SED parameters. This 
is a clear advantage of combining a range of data sets. Constraints 
on the SPIRE cirrus scale dependence and primary CMB and ACT 
radio source amplitudes are not shown as they are driven by the 
priors given in Section 3. 

The measured angular power spectra and best-fit model are 
shown in Figure I . The bandpowers shown include the corrections 
for source SEDs and bandpass filter transmission described in Sec- 
tion 3.6, as well as the best-fit calibration parameters. Figures 2 
and 3 show differential number counts from SPIRE and MIPS, re- 
spectively, with the binned counts likewise corrected for SED, filter 
and best-fit calibration. When comparing our model with other data 
sets, or predictions from other models, it is not sufficient to consider 
only the model curves plotted in Figures I to 3, which do not reflect 
the variation allowed by the calibration uncertainties. 

The remainder of this section deals, in turn, with the three 
parts of our model. We make comparisons with existing work and 
discuss a range of model extensions or modifications. 



4.2 Luminosity function 

While our constraint on the faint-end LF slope of qlf = 0.79tQ'42 
is consist ent with the value of ^ 1.2 typically obtained in the litera- 
ture (e.g. I Saunders et aLlll99(ll . BI 1), considerably lower values are 
also permitted. BlI found qlf = 1.223 ± 0.044; their constraint 
is tighter than ours only because of including very local constraints 
from IRAS and does not correspond to a tighter constraint on the 
behaviour of faint high-redshift sources. Given the mild preference 
shown for lower values of qlf, we consider allowing qlf to evolve 
as a power law in 1 + z, but with qlf fixed to 1 .2 at z = 0, and find 
that the data do indeed prefer a decrease in qlf with increasing z, 
but only at the 1.3(t level of significance. 

A range of values for (Jlf, which determines the bright-end 
behaviour, has been reported in the literature. Our constraint of 
0.34 ± 0.03 agree s fairly well with the 0.2010;^ obtained by 
ICaputi et al.l l l20Q7h from the 8 /im rest-frame luminosity function 
of star-forming galaxies al z — 1, for instance, but appears low 
compared to the 0.406 ±0.019 obtained by Bl I. We note, however, 
that (tlf is somewhat sensitive to the assumed degree of correlation 
between the stacked count uncertainties; if we neglect this corre- 
lation entirely, as in BIl, the constraint obtained is 0.38 ± 0.03, 
which agrees well with the B 1 1 result. As above, we consider al- 
lowing power-law redshift evolution of ctlf and find again a mild 
(1.2(t) preference for a decrease with increasing redshift. 

Including the deep GOODS-N counts from BI2 (which were 
not included in the baseline fit - see discussion in Section 3.9), as- 
suming a 50 per cent correlation between the errors on the GOODS 
counts in each redshift bin, as for the stacked COSMOS counts, 
yields tighter constraints of q:lf = 1-17 ± 0.12 and ctlf = 
0.29 ± 0.02. Being able to robustly extract deep number counts 
is clearly important for future faint-end LF constraints at high red- 
shift. 

Figure 4 shows the evolution of Lc with redshift, including la 
uncertainty contours. We also plot Lc{z) when additional freedom 
in the Lc and $c evolution is allowed (with a third order term in 
ln(l + z) for each - see equations 3 and 4). A preference for an in- 
crease in Lc beyond z ~ 2.5, and deviations from pure power-law 
evolution ml -\- z, are apparent in both cases. To further demon- 
strate the current data's high-redshift constraining power, we re- 
peated our fit with Lc {z) fixed to flatten into a plateau but otherwise 
be completely consistent with our best-fit model. Even allowing for 
higher-order terms in the <J>c evolution to compensate the lack of Lc 
freedom, this evolution is disfavoured at the 8, 6, and 3(t levels for 
plateaus at redshift 2.0, 2.5 and 3.0, respectively. 

It is possible that our Lc{z) results are being biased by our 
simplistic SED treatment, or failing to explicitly account for mul- 
tiple CIB source populations. Investigation of stacked count uncer- 
tainty correlations and different LF shape parametrizations is also 
clearly merited. We therefore do not attempt further interpretation, 
for instance inferring the star formation rate density, in the present 
work, however, our findings strongly suggest that joint fits to num- 
ber counts and power spectra can provide meaningful constraints 
on the evolution of the dust luminosity function and, ultimately, 
star formation rates at high redshift. 

Our approach is completely independent from existing anal- 
yses using luminosity function measurements of re solved sources 
(for recent examples u sing Spitzer data, see, e.g., iMagnelli et alj 
l20Ill ; |Patel et alj|2012l) . With future FIR, submm and microwave 
data improvements (see Section 5.7), we may expect analyses such 
as ours to yield constraints on high-redshift dust-obscured star for- 
mation that are both tighter and, thanks to directly probing the peak 
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Table 2. Marginalised parameter constraints 



Parameter Marginalised Description Equation 

constraint 

olp 0-79^q'42 power-law index describing faint-end slope of LF ("J" oc L^~" for low L) 2 

(TLF 0.34 ±0.03 LF spread 2 

log Lo I Lq 9.8 it 0.3 characteristic bolometric dust luminosity at redshift zero 3 

ei^ 5.2 ± 0.5 first order characteristic luminosity redshift evolution 3 

(^L ~0-98io'26 second order characteristic luminosity redshift evolution 3 

log $0 / g^il dex^^ Mpc^'* —1.9 it 0.2 LF normalisation at redshift zero 4 

— 2.5 ± 0.8 first order LF normalisation redshift evolution 4 

C,^ —0.3 it 0.5 second order LF normalisation redshift evolution 4 

To / K 24.1 it 1.4 temperature of dust component 'a' at z = 7 

zt 1.27^Q }g redshift for turn-on of component 'a' temperature evolution 7 

tT 0.75 it 0.10 component 'a' temperature evolution index 7 

Tf, /K 52 it 7 temperature of dust component 'b' (does not vary with redshift) 8 

P 1.78 it 0.11 dust spectral emissivity index 8 

/(, 0.46 it 0.10 contiibution of dust component '6' to total luminosity 8 

bo 0.84^Q :[g large-scale dusty source bias at redshift zero 10 

ei^ijjg 1.43 it 0.21 bias redshift evolution parameter 10 

Abias 1.61 it 0.19 source bias scale-dependence parameter 10 

fcc/Mpc~^ 4.9^Q'g small-scale bias truncation scale 10 

Acirr / Jy^ sr^i (1.3 it 1.0) X 10^ Galactic cirrus power at ^ = 2000 in SPIRE 1200 GHz spectrum 

^kSZ / CMB 5.3^2 4 amplitude of the kinematic Sunyaev Zel'dovich angular power spectrum aX £ = 3000 



Table 3. GIB source parameter coirelation matrix 
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or tail of the thermal dust emission, more robust, than those probing 
the rest-frame MIR. 



4.3 Spectral energy distribution 

We have been able to place tight constraints on our adopted ther- 
mal dust SED parameters, with two main results. Firstly, evolu- 
tion of the temperature of the 'a' dust component, which follows 

Ta = To (^Y+^y^ for z > zt, with To = 24.3 ± 1.4 K, 



ZT = 1.27l^;i5, and et = 0.75 ± 0.10, is strongly required 
(indicated by non-zero er). Secondly, an additional dust compo- 
nent, with redshift-independent temperature Tt = 52 ± 7 K, is also 
necessary - as shown by the significant preference for non-zero 
fb = 0.46 ± 0.10 (see equation 8). 

It seems possible that the evolution in Ta is due not to an evo- 
lution in the actual dust population, but to the increase in Lc{z); as 
discussed in Section 2.2, a positive correlation between temperature 
and luminosity has been observed in various studies. The apparent 



© 2012 RAS, MNRAS 000,[T]-?? 



Thermal dust emission in distant galaxies 13 



CD 
ft 



O 
ft 



a 

< 



10° 



10= 



10* 



10"= 



10= 



10* 









1 1 ' III. 

SPIRE 1200 GHz 












0^ 




10^ 

1 


10* 

1 1 1 1 



10 



10" 



10= - 



10* 



lo-" 



10^ 



10° 



10= r 



10* 



10^ 



10^ 




SPIRE 857 GHz 



10* 









1 1 1 


SPIRE 600 GHz 


I 


I 











lo-" 



10* 



Planck 857 GHz ■ 


s 


0^ 


10^ 






lv^anck217 GHz ; 






I \ 


< 

1 


> 



10'' 



lo-" 




Clustered power 
Poisson power 
Total power 



Multipole moment 



Figure 1. Angu lai' power s pectra of unresolved CIB sources from //erac/ie/-SPIRE jAmblard et ai]|201 ih . Planck jPlanck Collaboration et al]|201 ih . ACT 
jPas et al.ll201 ih , and SPT iReichardt et"ai]l2012h , with our best-fit model overplotted. The power spectra consist of a scale-independent Poisson term, and a 
clustered component that falls roughly as a power law. The bandpower uncertainties do not include uncertainties in map calibration (see Section 4.5). The best- 
fit Galactic cirrus power in the SPIRE spectra is approximately zero and not shown. SPIRE bandpowers from I < 2000 were not included in the fit (Section 3), 
but are shown here for comparison. The Planck bandpower uncertainties plotted include only the nominal beam uncertainties; additional freedom was allowed 
in the basehne model, as described in Section 3.7. The kSZ effect. Galactic cirrus and unresolved radio sources contribute to the ACT and SPT spectra but are 
highly subdominant to the dusty CIB source contribution. The best-fit primary lensed CMB has been subtracted from the ACT and SPT bandpowers (Section 
3.7.3). 



evolution in dust temperature may also be associated with evolu- 
tion in dust opa city. A more gener al form of the graybody equation 
is given by (e.g. lBlain et al ]l2003l) 

oc [1 - eiqi{-{v/v^f)\B^{T), (23) 

which asymptotes to the optically thin limit adopted in the base- 



line model for small v/uc. There is a strong degeneracy between 
r, /3 and Vc, and we find that the effect of introducing Uc may be 
largely reproduced by change s in Tg and /3 , at least for ~ 3000 
GHz (adopting the value from lDrainel2006l) . We may expect a sim- 
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Figure 2. Far-infrared differential number counts from Herschel-SPIRE iBethermin et al.l2012jl . obtained from a stacking analysis of the GOODS-N (squares) 
and COSMOS (triangles) fields, and by counting resolved sources in COSMOS (diamonds), with our best-fit model overplotted. The best-fit calibration values 
have been applied to the measured counts, as have corrections for source SED and the SPIRE bandpass filter transmission, as described in Section 3. The 
GOODS counts were not included in the baseline model (Section 3.9). 



ilar degeneracy to exist between parameters describing any redshift 
evolution in i/c and Ta- 

Regardless of the physical origin, the fact that the data require 
a non-zero er at a significance level of over 7a strongly suggests 
that some kind of SED evolution is required. We find that for the 
baseline parametrization, the evolution in temperature cannot be 
replaced with additional evolution terms in Lc or $c, or by intro- 
ducing evolution in emissivity index, /3, or allowing a different /3 
for the 'a' and components. 

We find that fitting for a power-law modification to the Wien 
tail, with free index, omir, as an alternative to the second dust 
component, yields similar results to the baseline model, provided 
evolution in the dust temperature is still allowed. The MIR SED 
index is constrained to be aum = 2.50 ± 0.23, somewhat higher 
than the value of 2 . adopted in recent work (e.g., [Hall et al. 2010; 
IShang et al.l |2012| : IViero et alj l2012bh . With this modified MIR 
model, the temperature of the dust at low redshift is higher than 
To in the baseline model by ~ 1.5a, and the emissivity index, /3, 
is ~ 1.5a lower Strikingly, the evolution parameters for the dust 
temperature are both in excellent agreement with those of the 'a' 
component in the baseline model, with temperature evolution again 
strongly required. 

The baseline and modified MIR SED models yield similar 
goodness-of-fit, although the baseline model is mildly favoured. 



Stronger discrimination between these models is not possible with 
the current data because of their similarity over the range of fre- 
quency probed (illustrated in Figure 5). Data from shorter wave- 
lengths, and an expanded SED model (including contributions 
other than thermal dust emission), will provide tighter constraints 
on the behaviour of the dust emission in the MIR range. We find no 
evidence for evolution in OfMiR when we adopt the modified MIR 
power law SED, and we likewise find no evidence for evolution in 
Tb in the baseline model. Allowing a modified MIR index in ad- 
dition to explicitly parametrizing a second temperature component 
does not lead to improvements in the model likelihood. 

While, as stated in Section 2, and discussed here, the extent to 
which our SED parameters can be associated with physical quan- 
tities is unclear, what our results demonstrate is that the counts 
and power spectra are capable of constraining the mean thermal 
dust SED at high redshift. It is also encouraging that our dust tem- 
perature and emissivity constraints show good general agreement 
with the r ~ 30 - 60 K and /3 ~ 1.5-2 obtai ned in SED fits 
to FIR-, submm- and mm-selected galaxies (e.g. 'Dunne & Eale j 
12001; Chapman et a l. 2005. ; .Coppin et al.,2006. : .Chapin et alj20Ill : 
loreve et alj20I2h . 

The SED constraints obtained in our model hinge on the use of 
data covering a wide range of frequency; in this work, the submil- 
limetre and microwave-band power spectra provide considerable 
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Figure 3. Far-infrare d differential number counts from Spi'fter-MIPS 
iBetliermin et alJlOlOl) . witli best-fit model overplotted. The counts are cor- 
rected for best-fit calibration, source SED and bandpass filter, as in Figure 
2 (see Section 3). 



improvement over the SPIRE and MIPS counts alone (see Sec- 
tion 5.5). Deep counts at lower frequencies than SPIRE (for in- 
stance, from the Submillimetre Common-User Bolometer Array-2 
- SCUBA-2 - Cosmology Legacy Surve}^) could likely fulfil this 
role to some extent as well, provided the effect of strong lensing on 
the bright-end counts can be robustly treated. 



4.4 Clustering properties 

The data strongly require scale-dependent dusty source biasing rel- 
ative to the linear matter power spectrum, parametrized through 
Abias = 1.61 ± 0.19 and K = 4.8lo 8 Mpc"\ A preference 
for scale-dependent bias remains at high significance even when 
the non-linear matter power spectrum (calculated using HALOFIT 
from the CAMB distribution) is used in our clustered power calcu- 
lation. 

Figure 6 shows the impact of this scale dependence as a func- 
tion of both I and fc; by ^ ~ 750, the scale dependence results 
in a factor of two difference compared to the linear matter power 
spectrum shape. The fact that the scale dependence has a significant 
impact even on these angular scale s - where the one-ha l o term of 
the halo model is subdom inant (e.g. lAmblard et alj20I ihlxia et alj 
l20I2l : IVieroetalJl20I2bh - suggests that it may not be valid to ne- 
glect the scale dependence of halo bias when calculating the two- 
halo halo model term, as discussed in Section I. 

Figure 7 shows our constraint on the evolution of the large- 
scale CIB source bias. Also shown is the bias of dark matter 



http://www.jach.hawaii.edu/JCMT/surveys/Cosmology.html 



10' 



10' 



10' 



baseline model 
extra evol. params y 




Plateau models: 

2 = 3.0; disfav. at 3cr 
2 = 2.5; disfav. at 6cr 
2 = 2.0; disfav. at 8cr 



_i_ 



1 



Figure 4. Redshift evolution of characteristic thermal dust luminosity, Lc- 
Model predictions are shown for the baseline model and when extra evolu- 
tion parameters are added to and <I>c. Three models that are consistent 
with the baseline predictions at low redshift, but feature a plateau in Lc 
at redshift 2.0, 2.5, and 3.0, are also shown. The ability of the counts and 
spectra to constrain high redshift evolution is demonstrated by the fact that 
these models are disfavoured at the 8, 6 and 3(T levels. 



haloes at several fixed masses as a function of redsh ift, using the 

[parametric fit to A^'-body simulations presented by [Tinker et alj 
2O10l) . Our baseline results are consistent with dusty sources in- 
habiting haloes of mass Afhaio ~ 10^^ Mq, which is higher than 
the values of ~ lO'^ Mp) associated with optimal star formation 
from recent studies (e .g. IWang et al]|20I2l : iBethermin et al ]|2012bl : 
iBehroozi et aLll2012b||3) . Large number of fainter sources occupy- 
ing massive groups and clusters may lead to an increase in the ef- 
fective bias we have measured. We also find that lower values of bo, 
corresponding to a lower characteristic halo mass, are allowed if we 
introduce a correlation between CIB source bias and dust luminos- 
ity (see Section 2.4). As an example, we consider a correlation of 
the form 



(fogal(fc,Z;)>(2))t:ut = (6gal(fc,2)) 



(>(2))cut 



+ 7corr(l + Z)X 



Sou 



dL ^{L,z) ( S, 



dS^ Lin 10 V Sc 



5. 



(24) 



where 7corr is a free parameter, which may be positive or negative, 
and Sc,v ~ Sv{Lc), such that the term in parenthesis increases as 
the relative importance of the more luminous sources increases. We 
find a mild preference for a positive value of 7corr, which leads to 
a lower inferred bias, shown in Figure 7. As stated in Section 2.4, 
the data show only a mild (2a) preference for a &gai — L correlation 
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Figure 5. Comparison of the best-fit dust SED from the baseline model and 
an alternative model featuring a power-law modification to the MIR Wien 
tail instead of a second dust component, both shown at z < , where the 
SED is not redshift-dependent, for the same integrated thermal dust lumi- 
nosity. The shapes of the SEDs from the two models are extremely similar 
over a wide range of frequency, and the data considered do not strongly 
prefer one over the other. 
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Figure 6. Dusty galaxy power spectrum, Pgai, as a function of multipole 
moment and wavenumber (related by fc = + l/2)/x for small angles). 
We show the linear matter power spectra, scaled to match the large-scale 
Pgal behaviour, to illustrate the impact of the scale-dependent source bias, 
which is important for fc > 0.1 Mpc^^. At z ~ 1, where the contribution 
to the clustered power peaks, the scale dependence results in a doubling of 
fgai by t ~ 750. 
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Figure 7. Evolution of large-scale CIB source bias: we show Icr contours 
from both our basehne model, and with the addition of the b — L correlation 
from equation 24. A positive b — L correlation is mildly preferred, leading 
to lower bias values. The bias of dai'k matter haloes of several fixed masses 
ar e also shown as a fu nction of redshift, using the fitting functions provided 
bv lTinkeret"ai]i2010l) . 



of this form in terms of goodness-of-fit, and the introduction of the 
correlation does not impact significantly on LF or SED constraints. 

Further work, involving more detailed modelling of the CIB 
source clustering, perhaps within a halo model framework, is re- 
quired to explore the relation between dust luminosity and halo 
mass in more detail. 



4.5 Calibration 

Constraints on the photometric calibration parameters, which were 
included as additional free parameters in our fitting (Section 3.5) 
are shown in Table 4. The fact that the marginalised uncertainties 
from our fit are smaller than the nominal uncertainties is a con- 
sequence of fitting to a large number of data sets with substantial 
frequency overlap. As discussed in Section 3, what we are calling 
calibration parameters may also be absorbing uncertainties relating 
to bandpass filter transmission profiles and the filter response to the 
source SEDs. We therefore do not expect the /cai constraints to be 
consistent with unity in every case. 

Table 4 shows the constraints on the additional SPIRE beam 
area correction factors, /beam, introduced as described in Section 
3.7. Also shown are the constraints obtained when neither these fac- 
tors, nor the additional Planck beam uncertainty, are allowed ("No 
beam corrections" column). Removing the additional beam-related 
freedom does not significantly affect the calibration values, except 
that the SPIRE spectrum calibrations absorb the effect of the beam 
area correction and lie further from unity. The marginalised CIB 
parameter constraints also change by less than 0.5a in all cases. 
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Table 4. Marginalised constraints on calibration parameters 

Data set Baseline No beam Nominal 

model corrections 




10^ 

Frequency [GHz] 



Figure 8. Mean CIB intensity spectrum, measured using FIRAS 
(Lagache et al. 199^ gr een line) and from an analysis of IRAS and MIPS 
data at 100 and 160 /^m jPenin et al .l2012bl . red diamonds). The solid black 
line shows our mean model prediction, and the dashed lines the Icr confi- 
dence contours. 

surements of the integrated CIB SED from FIRAS l lLagache et alj 
We also show measurements of the mean CIB inte nsity at 
160 and 100 ^im from IRAS and MIPS l lPenin et alj|2012d) . 

Our model is in fairly good agreement with the CIB intensity 
measurements, with an underestimation of ~ Icr around the peak of 
the CIB emission. In principle, we could have used the FIRAS data 
to constrain the model, along with the spectra and number counts; 
we did not consider this because the Planck 857 and 545 GHz maps 
were calibrated from the FIRAS measurements, and thus any sys- 
tematic present in the FIRAS data could have effectively entered 
our modelling twice. 



5.2 Power spectra at lower frequencies 

As stated in Section 3, we limited our analysis to spectra contain- 
ing negligible contribution from the Sunyaev Zel'dovich effect or 
its cross-correlation with CIB sources. A good test of our model as- 
sumptions is to extrapolate the power spectrum predictions to lower 
frequencies. We compare our predicti ons with the most rece nt CIB 
power spectrum constraints from SPT l lReichardt et al .120121) in Ta- 
ble 5. We compare to the SPT measurements obtained when a cor- 
relation between the t SZ effect and CIB sources is allowed, since 
lAddison et al. found that such a congelation may exist at 

the 10-20 per cent level in units of power. 

There is fairly good consistency between our predictions and 
the SPT measurements, although our predictions for Poisson power 
are somewhat higher than the SPT results. Our constraints, at least 
on the clustered power amplitude, are competitive with the direct 
measurements. Extending our model to include lower frequency 
{v < 220 GHz) ACT, SPT, and Planck data may therefore im- 
prove constraints on both our dust emission model, and secondary 
anisotropics like the Sunyaev Zel'dovich effect. 
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t an additional 2 per cent has been added to the Planck calibration 
uncertainties (Section 3.6.1) 
t SPT calibration is in units of power, otherwise units of flux or CMB 
temperature 

* we enforce congelation between these calibration values, with 
band-to-band covariance of 0.05^ (Section 3.5) 

We allowed additional freedom in the SPIRE and Planck 
be am treatment in the baseline mo del, motivated by discus sions 
in lPlanck Collaboration et alj ( 1201 ih and lViero et~ai] ( l2012bl) . but 
have demonstrated that this is not significantly affecting or bias- 
ing our results. We further note that removing either the SPIRE 
or Planck spectra from the fit entirely does not qualitatively mod- 
ify the shape of the luminosity function evolution or the require- 
ment for SED evolution and the second temperature component 
discussed earlier in this section. Additionally, the preference for 
the low Planck 545 GHz calibration remains even when the SPIRE 
spectra are not included. 



5 DISCUSSION 

In this section we compare our model's predictions with several 
measurements outside those used to constrain it, and discuss some 
of its implications for future work. 

5.1 Integrated CIB intensity 

The integrated CIB intensity, clustered anisotropy power and Pois- 
son anisotropy power each receive different contributions from dif- 
ferent populations of sources, with the brightest sources making a 
larger relative con tribution to the Poi sson p ower, for instance (se e 
discussions in e.g. lHaiian et alj|2012l . A12. lAddison et alj[2oT2bh . 
Our model correctly reproduces the anisotropy power spectra over 
a range of frequency; an obvious question is whether it also repro- 
duces the integrated CIB intensity. We show our mean model pre- 
dictions, with Ict uncertainty contours, in Figure 8, along with mea- 
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Table 5. Power spectmm predictions for lower frequencies (units of I 
l)C,/27r|,=3ooo,m /^K^CMB) 

Prediction Measured^ 



150 GHz 



95 GHz 



Poisson 
Clustered 
Poisson 
Clustered 



9.6 ±0.9 
7.0 ±0.8 
1.37 ±0.21 
0.97 ±0.17 



8.1 ±0.5 
6.7 ±0.7 
1.04 ±0.15 
0.87 ±0.17 



t SPT constraints from lReichardt et al] j2012l) with tSZ-ClB coiTelation 
allowed 



5.3 Kinematic Sunyaev Zel'dovich constraints 



c q+2.2 



The kSZ constraint obtained in our baseline fit, Aksz 
/^K^, represents a preference for a non-zero kSZ contribution at 
slightly more than the la significance level, which remains even 
if no upper limit is imposed on the kSZ ampl itude. The ampli- 
tude is co nsistent with recent predictions (e.g. IShaw et al.l [20121 ; 
iMesinger e t al. 2012), however the current uncertainty is too large 
for our result to be useful for constraining kSZ models. We find 
that the kSZ amplitude is not highly correlated with any of the dust 
model parameters, meaning that there is no single quantity which, 
if better known, would lead to a substantial improvement in the kSZ 
constraint. While constraining the kSZ power spectrum from data 
free from the tSZ is an attractive idea, it would appear that, at least 
currently, this is not feasible, given the small size of the kSZ power 
relative to the dust at 220 GHz. 
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Figure 9. Two-dimensional parameter constraints with contours enclosing 
68 and 95 per cent of the MCMC samples. Constraints on SED parame- 
ters (such as the low-redshift temperature of the 'a' component. To, and 
the parameter describing the temperature evolution, tT - see equation 7) 
are considerably improved by fitting the counts and clustered power over 
fitting the SPIRE and MIPS counts alone. Requiring the model to also fit 
the Poisson power spectrum further improves the SED constraints, and also 
significantly improves constraints on the large-scale bias of the CIB sources 
(given at 2: = by fco ) and the redshift evolution of the bias, parametrized 
by ebias- 



5.4 The signal in the shot-noise 

An interesting result from our work is the constraining power that 
is carried by the Poisson component of the high-resolution spectra 
(mainly SPIRE, but also ACT and SPT). To illustrate this. Figure 
9 shows two-dimensional parameter constraints, with contours en- 
closing 68 and 95 per cent of the MCMC samples, for several model 
parameters. We show constraints from three models: the baseline 
model, a model constrained from the SPIRE and MIPS number 
counts alone, and a model constrained by the counts plus clustering 
power alone. In this third model, the Poisson power in each of the 
nine power spectra is allowed to vary as an additional free param- 
eter. The results of this comparison can be summarised as follows: 
including the clustered power leads to a significant improvement in 
SED parameter constraints over the counts alone, and inclusion of 
the Poisson power in the fit leads to moderate further improvements 
in the SED parameters constraints, and large improvements in the 
bias parameter constraints, in particular dramatically reducing the 
correlation between 6o and ebias- 

We also find that the model likelihood is not sufficiently im- 
proved by treating the Poisson levels independently to warrant a 
preference for this approach on goodness-of-fit grounds (Ax^ = 
15 for 9 fewer degrees of freedom, corresponding to a I.3ct prefer- 
ence). Furthermore, there is good consistency between the results 
from the three models for every parameter. This is a useful test of 
at least two features of our model. Firstly, it suggests that our re- 
sults are not significantly biased by our treatment of bright source 
masking prior to power spectrum calculation (see Section 2.5). Sec- 
ondly, it suggests that the fact that we have not explicitly accounted 
for a dusty galaxy duty cycle in our model (which would have most 
impact on the pre dicted Poisson power - see, e.g., discussion in 
IShangetal.ll2012h is not significantly affecting our results. 



5.5 Redshift distribution of power and intensity 

Figure 10 shows the redshift distribution of clustered and Poisson 
anisotropy power, and CIB intensity, as a function of frequency, 
for our best-fit model parameters. Different levels of bright source 
removal is responsible for the behaviour of the Poisson power at 
low redshift, with notably fewer bright sources bring masked in the 
Planck compared to SPIRE spectra. Our model predicts that, as ex- 
pected, lower freque ncies receive greate r relative contribution from 
higher redshifts (e.g. iKnox et al]|200ll) . but, interestingly, there is 
only moderate evolution in the distribution over the submm and mi- 
crowave bands. Physically, this corresponds to a high degree of co- 
herence in the dust emission across this frequency range. Defining 
the degree of cross-correlation between frequency v\ and V2 as the 
ratio of the cross-spectrum power to the square root of the product 
of the auto-spectra, Ci^n^^^j Ci^i,^ Ci,v2, we predict a correlation 
in units of power of around 90 per cent between the Planck 857 and 
217 GHz bands, for example. 

This high degree of correlation would suggest that much of the 
dust that acts as a foreground contaminant for current CMB temper- 
ature cosmology (hampering detection of the kSZ power spectrum, 
for instance) may be removed either by direct cross-correlation of 
microwave maps with maps at higher frequencies, or understood 
through a joint fit with the high-frequency, dust-dominated data. 

The most direct test of this prediction is comparison with mea- 
sured power spectra from cross-correlating microwave and submm 
maps. We repeated our model fitting including cross-spectra pre- 
sented by iHaiian et alj ( |2012[) . calculated by correlating BLAST 
maps at 250, 350 and 500 fim with the ACT 218 GHz map, and 
found that our baseline model provides a good fit (Ax^ ~ 25 for 
27 additional bandpowers, treating the BLAST map calibration pa- 
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rameters as described in Section 3.5). Ttiis is consistent with tlie 
iiigh degree of colierence pref erred by most of the ACT / BLAST 
cross-spectra in the analyses o f lHaiian etal] ( l2012h and A 12, how- 
ever, the BLAST / ACT cross-spectra bandpowers are noisy, and 
future cross-correlations (for instance from Planck, or from cross- 
correlating SPIRE maps with ACT or SPT) will provide a stronger 
test of our prediction. 

Since the evolution in dCi/dz with changing frequency is 
more significant for the Poisson than the clustered power, we fur- 
ther predict that the degree of correlation between different bands 
is ^-dependent, with the correlation highest on large angular scales, 
where the Poisson is subdominant. Our predictions for the scale 
dependence of the SPIRE cross-c orrelations a r e in goo d agreement 
with the recent measurements of Iviero et al.l ( 1201 2bl) . We predict 
cross-correlations of 0.92-0.96, 0.85-0.91 and 0.95-0.98 for the 
SPIRE 1200x860, 1200 x 600 and 860 x 600 GHz spectra, respec- 
tively, w ith the range indicat ing the variation from £ ~ 2 x 10'* to 
I ~ 10^. IVieroet"^ ( l2012hl) find cross-correlations of 0.95±0.04, 
0.86 ± 0.04, and 0.95 ± 0.03 for the same three cross-spectra (also 
see their Figure 10, which shows the trend of decreasing correlation 
with increasing I). 



5.6 Model limitations and future work 

A large assumption in our modelling is that the dusty sources can 
essentially be treated as a homogenous population, whose aver- 
age properties vary smoothly over time. Much of our constraining 
power is in the form of integrated quantities - the power spectra 
- meaning our results are probably not especially sensitive to rare, 
extreme objects (which could distort results obtained from a flux- 
limited study), but it is not clear that our model would correctly 
capture a sharp transition in either dust or clustering properties, 
which feature s in some existing work. In the model develope d by 
iGranato et all ( l2004h . lLaDi et al.l ( 1201 ih . and lXia et al] ( |2012|) . for 
instance, highly-clustered, massive, proto- spheroidal galaxies dom- 
inate the dusty galaxy population at high redshift, but cease to exist 
at 2 < 1, where merger-driven starbursts and (eventually) spiral 
galaxies take over. Allowing more explicitly for multiple source 
populations, with separate emission or clustering properties, is a 
clear improvement we can consider in future work. 

Some possible limitations of our SED parametrization have 
already been discussed (Sections 2.2 and 4.3). The Atacama Large 
Millimeter/submillimeter Array (ALM/13) will greatly improve our 
understanding of the processes relating to dust-enshrouded star for- 
mation on galaxy scales, and may enable us move on from the gray- 
body approximation to a more physical dust emission treatment. 
Additionally, extending our model to include higher frequencies, 
incorporating non-thermal dust contributions to the SED and bolo- 
metric IR luminosity, will facilitate comparison with infrared lumi- 
nosity density and star formation rate constraints from other work. 

The first release of Planck maps, covering the entire sky 
(rather than just the ~140 deg'^ used for the spectra we have in- 
cluded here), is expected in early 2013. Spectra from these maps, 
in add ition to futur e Herschel, ACT, SPT , ACTPol ([Niemack et aL 
I2OIOI) and SPTpol dMcMahon et al.ll2009l : lAustermann et arj2012h 
spectra, will significantly improve constraints over those con- 
sidered here. Future number counts, from additional SPIRE 
fields, as well as, for example, SCUBA-2, will also give much- 
needed improved constraints on the properties of fainter high- 
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redshift sources. The angular correlation function of flux-limited 
samp l es of CIB sourc es (e.g. ICoorav et alj [ioiol ; iMaddox et alj 
l201(]|:lGuo etalJl201ll) . from, for instance, the //eric/je/-ATLAS 
( lEalesetalJl2010l)~ or SCUBA-2 will, when used in conjunction 
with the angular power spectrum, allow more detailed investigation 
of the 6gai — L congelation. 

A host of CIB cross-correlations, using other tracers of large- 
scale structure, such as reconstructed CMB or galaxy weak lensing 
deflection maps, or optical data from the SDSS, are expected in 
the coming months. Being sensitive to different combinations of 
model parameters from the counts and CIB spectra, these statistics 
will be invaluable for furthering out understanding of extragalactic 
dust emission. The framework presented in Section 2 can be easily 
used to either predict these correlations, or extended to use them to 
improve constraints. Our model is particularly well-suited to pre- 
dicting or interpreting the CMB lensing or integrated Sachs-Wolfe 
correlations with the CIB, where much of the measured signal will 
be on large angular scales, and relatively insensitive to the details 
of the small-scale clustering. 

The primary motivation for studying extragalactic dust emis- 
sion is to ultimately constrain the star formation history. Detec- 
tors onboard, for example, Spitzer and WISE, probe re-processed 
starlight in the near- and mid-infrared, complementing the longer- 
wavelength data we have used. T he near- and far-inf rared emission 
are known to be correlated (e.g. iBond et al]|2012l) . and incorpo- 
rating data at shorter wavelengths in our model will both improve 
understanding of this correlation, and lead to a more complete view 
of obscured star formation. 



6 CONCLUSIONS 

We have simultaneously constrained the luminosity function, spec- 
tral energy distribution and clustering properties of dusty CIB 
galaxies using a combined fit to number counts and angular power 
spectra from five current instruments - Spitzer, Herschel, Planck, 
ACT and SPT. Our analysis is based on fairly simple, phenomeno- 
logical, parametrizations, but has yielded a number of important 
results: 

(i) evolution in the dust SED, parametrized here through evolu- 
tion in graybody dust temperature, appears to be required, 

(ii) the data show a strong preference for either a second gray- 
body component or modified MIR tail behaviour, 

(iii) joint fits to FIR, submm and microwave counts and spectra 
can constrain the evolution of the dusty source luminosity function 
at high redshift (2 > 2.5 - 3), 

(iv) pure power-law CIB source clustering is moderately dis- 
favoured by combining current Planck and SPIRE spectra, 

(v) we predict a high degree of submm / microwave band CIB 
correlation (about 90 per cent in units of power between Planck's 
217 and 857 GHz bands, for example), and 

(vi) the Poisson, shot-noise, component of the power spectrum 
can be used to improve CIB source constraints. 

None of the data sets considered in this work are individually 
capable of yielding constraints competitive with those from the 
combined fit. The basic approach we put forward - combining 
data from a range of instruments, with a rigorous treatment of data 
and modelling uncertainties - will thus be of great importance for 
realising the potential of future data sets. 
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Figure 10. Redshift distribution of clustered and Poisson CIB anisotropy power, and CIB intensity, normalised such that the area under each curve is unity, for 
our best-fit model. The redshift distribution of the clustered power is shown at £ = 3000, although the redshift distribution does not vary strongly with scale 
over the range of scales where the clustered power contributes significantly to the total spectrum. A few per cent of the Poisson power (particularly for Planck, 
where fewer bright sources are removed) is predicted to originate from low-redshift sources that are below the source removal limit, although in each case the 
power remains convergent. The high degree of overlap between the curves at different frequency corresponds to a high degree of coherence in the dusty galaxy 
emission over the submillimeter and microwave bands. This is an important prediction of our model. 
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